NB : Ce script est un document de travail.
Il cherche à objectiver, par l’utilisation de méthodes d’analyse spatiale, l’« effet de grappe » repéré par la littérature (Bideau 2019)
Il est mis à disposition dans une logique de science ouverte.
Ce travail s’inscrit dans le cadre d’une étude plus générale sur les communes nouvelles (et plus particulièrement d’un addendum au chapitre IV de la thèse) :
https://cv.hal.science/gabriel-bideau
Licence CC-BY-NC-SA.
Il est possible d’accéder au code de ce Markdown ici : https://gbideau.github.io/CN_effet_de_grappe/effet_de_grappe.Rmd
Les données utilisées pour jouer le code sont regroupées ici : https://gbideau.github.io/CN_data/
Ne pas hésiter à contacter l’auteur (gabriel.bideau@gmail.com) pour toute question.
# Librairies utilisées
library(sf)
library(cartography)
library(mapsf)
library(corrplot)
library(cowplot)
library(MTA)
library(readxl)
library(ggplot2)
library(FactoMineR)
library(factoextra)
library(cluster)
library(reshape)
library(reshape2)
library(flows)
# NB : Pour le package flows, la version la plus récente est disponible ici :
# remotes::install_github("rCarto/flows") # ou # install.packages("mapsf")
# Pour obtenir une version plus ancienne (celle utilisée ici) : https://cran.r-project.org/src/contrib/Archive/flows/
# install.packages("packages/flows_1.1.1.tar.gz", repos=NULL, type="source")
library(sp)
library(knitr)
library(condformat) # https://cran.r-project.org/web/packages/condformat/vignettes/introduction.html
library(units)
library(stringr)
# library(dplyr)
library(questionr)
library(spdep) # Pour les matrices de contiguïté
library(rgeoda) # Pour les matrices de contiguïté
# Liste pour installer les packages si besoin :
# sf cartography mapsf readxl foreign dplyr flextable knitr stringr units condformat forcats ggplot2 rstatix questionr corrplot gtsummary broom GGally effects forestmodel ggeffects labelled cowplot spdep rgeoda
load("data/refdata.Rdata")
geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE)
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE)
Il s’agit ici de creuser l’« effet de grappe », ce qu’on appelle parfois l’« effet consultant ».
À l’aide de Bellefon, Loonis, and Le Gleut (2018).
library(spdep)
test_geom <- merge(geom2011, df2011[, c("CODGEO", "CODE_DEPT", "COM_NOUV", "P09_POP", "LIBGEO")], by = "CODGEO")
test_geom <- subset (test_geom, CODE_DEPT == "49" )
test_queen <- poly2nb(test_geom, queen = TRUE)
test_rook <- poly2nb(test_geom, queen = FALSE)
# On change de type d'objet pour faciliter les représentations
test_geom <- as(test_geom, "Spatial")
# Représentation graphique des deux manières de calculer les voisinages.
plot(test_geom, border="lightgray")
plot(test_queen, coordinates(test_geom),add=TRUE,col="red")
plot(test_rook, coordinates(test_geom),add=TRUE,col="blue")
# Réalisation d'une liste de poids
test_liste <- nb2listw(test_queen, zero.policy = TRUE)
# NB : zero.policy :default NULL, use global option value; if FALSE stop with error for any empty neighbour sets, if TRUE permit the weights list to be formed with zero-length weights vectors https://r-spatial.github.io/spdep/reference/nb2listw.html
# Pour avoir une représentation graphique des communes nouvelles
test_geom2 <- merge(geom2011, df2011[, c("CODGEO", "CODE_DEPT", "COM_NOUV", "P09_POP", "LIBGEO")], by = "CODGEO")
test_geom2 <- subset (test_geom2, CODE_DEPT == "49" )
plot(subset(test_geom2[, c("CODE_DEPT", "COM_NOUV")], CODE_DEPT == "49"))
# Pour extraire les entités sans liens :
test_liste$neighbours
## Neighbour list object:
## Number of regions: 363
## Number of nonzero links: 2012
## Percentage nonzero weights: 1.526915
## Average number of links: 5.5427
isol <- test_geom2[c(5896, 7662, 10958, 10959, 10960, 11024, 21454, 21470, 21471, 21472, 21473, 33974),]
rm(isol, test_queen, test_rook)
À l’échelle du département du Maine-et-Loire, des différences notables en fonction du voisinage choisi.
Officiellement, le contact par un point autorise à fusionner en commune nouvelle donc choix du voisinage “queen”.
Cf. la réponse du ministère de l’intérieur à une question écrite au Sénat sur la question de la continuité territoriale en 2009 : https://www.senat.fr/questions/base/2009/qSEQ090609011.html.
Cf. Bellefon and Salima (2018) p. 56 sq. pour l’explication et 62 pour le code.
test_geom$COM_NOUV2 <- dplyr::if_else(test_geom$COM_NOUV == "OUI", 1, 0)
moran.plot(test_geom$COM_NOUV2, test_liste, labels=FALSE , xlab="variable",ylab="moyenne des voisins")
moran.plot(test_geom$P09_POP, test_liste, labels=FALSE , xlab="variable",ylab="moyenne des voisins")
Cf. Bellefon and Salima (2018) p. 62 sq. pour l’explication et 64 pour le code.
L’analyse des statistiques des join count observe, pour une variable binaire, s’il y a une autocorrélation spatiale de la présence d’individus partageant ou non une caractéristique. “on considère une variable binaire qui représente deux couleurs, Blanc (B) et Noir (N) de sorte qu’une liaison puisse être qualifiée de Blanc-Blanc, Noir-Noir ou Blanc-Noir. On observe :
— une autocorrélation spatiale positive si le nombre de liaisons Blanc-Noir est significativement inférieur à ce que l’on aurait obtenu à partir d’une répartition spatiale aléatoire ;
— une autocorrélation spatiale négative si le nombre de liaisons Blanc-Noir est significativement supérieur à ce que l’on aurait obtenu à partir d’une répartition spatiale aléatoire ;
— aucune autocorrélation spatiale si le nombre de liaisons Blanc-Noir est approximativement identique à ce que l’on aurait obtenu à partir d’une répartition spatiale aléatoire.” (p. 62-63, souligné dans le texte)
test_geom$COM_NOUV3 <- as.factor(test_geom$COM_NOUV)
class(test_geom$COM_NOUV3)
## [1] "factor"
joincount.test(test_geom$COM_NOUV3, test_liste, zero.policy=TRUE, alternative="greater",
sampling="nonfree", spChk=NULL, adjust.n=TRUE)
##
## Join count test under nonfree sampling
##
## data: test_geom$COM_NOUV3
## weights: test_liste
##
## Std. deviate for NON = 12.114, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 44.663005 26.113260 2.344626
##
##
## Join count test under nonfree sampling
##
## data: test_geom$COM_NOUV3
## weights: test_liste
##
## Std. deviate for OUI = 11.408, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 89.519769 69.613260 3.045139
Pour les communes nouvelles, on observe une autocorrélation positive très nette et significative. L’autocorrélation est d’ailleurs également significtive pour les communes non fusionnantes, même si la déviation est moins nette.
Mesures locales d’association spatiale : LISA (Local Indicators of Spatial Association).
Cf. Bellefon and Salima (2018) p. 65 sq. pour l’explication et 69 pour le code.
“Anselin (ANSELIN 1995) développe les notions introduites par Getis et Ord en définissant des indicateurs d’autocorrélation spatiale locale. Ceux-ci doivent permettre de mesurer l’intensité et la significativité de la dépendance locale entre la valeur d’une variable dans une unité spatiale et les valeurs de cette même variable dans les unités spatiales environnantes.” (p. 65)
Fonction ‘localmoran’ du package ‘spdep’. Ii : local moran statistic ; E.Ii : expectation of local moran statistic ; Var.Ii : variance of local moran statistic ; Z.Ii : standard deviate of local moran statistic ; Pr() : p-value of local moran statistic (https://www.rdocumentation.org/packages/spdep/versions/1.1-8/topics/localmoran).
Pour analyser correctement la représentativité des résultats, plusieurs méthodes proposent d’ajuster la p-value :
“On considère que cette méthode [celle de Bonferroni] ne donne de bons résultats que lorsque le nombre de tests réalisés est petit. […] La méthode d’ajustement de Holm conduit à un plus grand nombre de clusters significatifs que la méthode de Bonferroni. Elle lui est donc le plus souvent préférée. Cependant, cette méthode se concentre aussi sur la détection de l’existence d’au moins un cluster dans toute la zone. […] La méthode du False Discovery Rate (FDR) a été introduite par BENJAMINI et al. 1995. Avec cette méthode, le risque de juger - à tort - un cluster comme significatif est plus élevé, mais inversement le risque de juger - à tort - un cluster comme non significatif est plus faible.” (p. 68) “La méthode de Holm diminue en effet le risque de conclure à tort à l’existence d’une autocorrélation spatiale locale. En revanche, cette méthode augmente le risque de passer à côté d’un cluster local.” (p.69)
results_LISA<- as.data.frame(localmoran(test_geom$COM_NOUV2,test_liste,zero.policy=TRUE))
#Calcul des p-values ajustées
results_LISA$pvalue_ajuste_bonferroni<- p.adjust(results_LISA$`Pr(z != E(Ii))`,method="bonferroni")
results_LISA$pvalue_ajuste_holm<- p.adjust(results_LISA$`Pr(z != E(Ii))`,method="holm")
results_LISA$pvalue_ajuste_fdr<- p.adjust(results_LISA$`Pr(z != E(Ii))`,method="fdr")
test_geom3 <- cbind(test_geom2, results_LISA)
# Quelles communes ont une p-value élevée selon les trois méthodes ?
par(mar = c(0,0,0,0), mfrow = c(2,2))
# plot(st_geometry(dep), col = NA)
choroLayer(x = test_geom3 , var = "pvalue_ajuste_holm",
method = "quantile", nclass = 6,
# col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
col = carto.pal(pal1 = "blue.pal", n1 = 6),
border = NA,
legend.pos = "topleft", legend.values.rnd = 2,
legend.title.txt = "Pvalue ajustée : méthode de Holm")
choroLayer(x = test_geom3 , var = "pvalue_ajuste_bonferroni",
method = "quantile", nclass = 6,
# col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
col = carto.pal(pal1 = "blue.pal", n1 = 6),
border = NA,
legend.pos = "topleft", legend.values.rnd = 2,
legend.title.txt = "Pvalue ajustée : méthode de Bonferroni")
choroLayer(x = test_geom3 , var = "pvalue_ajuste_fdr",
method = "quantile", nclass = 6,
# col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
col = carto.pal(pal1 = "blue.pal", n1 = 6),
border = NA,
legend.pos = "topleft", legend.values.rnd = 2,
legend.title.txt = "Pvalue ajustée : méthode FDR")
typoLayer(x = test_geom3, var = "COM_NOUV",
col=c("red", "blue"), border = NA)
layoutLayer(title = "Titre", theme = "red.pal",
author = "G. Bideau, 2022.",
sources = "")
# On choisit la méthode FDR et on ne conserve que les valeurs qui nous paraissent significatives au vu de la p-value ainsi ajustée
test_geom3$Ii_retenu <- ifelse(test_geom3$pvalue_ajuste_fdr <= 0.1, yes = test_geom3$Ii, no = NA)
par(mar = c(0,0,0,0), mfrow = c(1,2))
# On cartographie le résultat de l'indice de Moran
choroLayer(x = test_geom3 , var = "Ii_retenu",
method = "quantile", nclass = 6,
# col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
col = carto.pal(pal1 = "blue.pal", n1 = 6),
legend.pos = "topleft", legend.values.rnd = 2,
legend.title.txt = "Indices de Moran significatifs",
legend.nodata = "p-value non significatif (>0,1 méthode FDR)")
typoLayer(x = test_geom3, var = "COM_NOUV",
col=c("red", "blue"), border = NA, legend.pos = "topleft")
par(mar = c(0,0,0,0), mfrow = c(1,1))
Présentation du package ‘rgeoda’ : https://geodacenter.github.io/rgeoda/articles/rgeoda_tutorial.html
Sur les matrices de contiguïté : https://geodacenter.github.io/workbook/4a_contig_weights/lab4a.html
Sur l’auto-corrélation : http://geodacenter.github.io/workbook/6a_local_auto/lab6a.html
library(rgeoda)
test_geom <- merge(geom2011, df2011[, c("CODGEO", "CODE_DEPT", "COM_NOUV", "P09_POP")], by = "CODGEO")
# test_geom <- subset (test_geom, CODE_DEPT == "49" )
test_geomCN <- merge(geom_new, df_new, by = "CODGEO_new")
test_geomCN <- subset(test_geomCN, COM_NOUV == "OUI")
queen_w <- queen_weights(test_geom, order=1, include_lower_order = FALSE, precision_threshold = 0)
summary(queen_w)
## name value
## 1 number of observations: 36208
## 2 is symmetric: TRUE
## 3 sparsity: 0.000164767983149941
## 4 # min neighbors: 0
## 5 # max neighbors: 29
## 6 # mean neighbors: 5.96591913389306
## 7 # median neighbors: 6
## 8 has isolates: TRUE
# Pour savoir si certaines entités sont séparées des autres
has_isolates(queen_w)
## [1] TRUE
NB : Sous-section qui peut poser problème lors d’un knit si n’est pas executée en même temps que les précédentes (cela ne fonctionne pas bien si les sous-sections précédentes sont en cache).
test_geom$COM_NOUV2 <- dplyr::if_else(test_geom$COM_NOUV == "OUI", 1, 0)
variable <- "COM_NOUV2"
df_variable <- test_geom[variable]
lisa <- local_moran(queen_w, df_variable)
# To get the values of the local Moran:
lms <- lisa_values(gda_lisa = lisa)
# To get the pseudo-p values of significance of local Moran computation:
pvals <- lisa_pvalues(lisa)
# To get the cluster indicators of local Moran computation:
cats <- lisa_clusters(lisa, cutoff = 0.05)
# Labels
lbls <- lisa_labels(lisa)
test_geom$LISA <- lms
test_geom$pvals <- pvals
test_geom$cats <- cats
# On rend les catégories plus lisibles en remplaçant le numéro par l'intitulé
test_geom$cats <- as.factor(test_geom$cats)
num_levels <- as.numeric(levels(test_geom$cats))
levels(test_geom$cats) <- lbls[num_levels + 1]
print(levels(test_geom$cats))
## [1] "Not significant" "High-High" "Low-High" "High-Low"
## [5] "Isolated"
test_geom$LISA_retenu <- ifelse(test_geom$pvals <= 0.1, yes = test_geom$LISA, no = NA)
par(mar = c(0,0,0,0), mfrow = c(1,2))
couleurs <- c("#f7f7f7", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33")
# Cf. https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6 pour composer les palettes
typoLayer(x = test_geom , var = "cats",
# col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
col = couleurs[1:length(levels(test_geom$cats))],
border = NA,
legend.pos = "topleft",
legend.values.order = levels(test_geom$cats),
legend.title.txt = "Indices de Moran significatifs\nsur la variable de « commune nouvelle »\n(pvalue < 0,05)")
typoLayer(x = test_geom, var = "COM_NOUV",
legend.values.order = c("OUI", "NON"),
col=c("red","blue"),
border = NA,
legend.pos = "topleft",
legend.title.txt = "Communes fusionnantes")
# plot(test_geomCN$geometry, col = NA, border = "black", lwd = 1, add = TRUE)
par(mar = c(0,0,0,0), mfrow = c(1,1))
table(test_geom$cats, test_geom$COM_NOUV)
##
## NON OUI
## Not significant 31245 992
## High-High 0 1576
## Low-High 2380 0
## High-Low 0 3
## Isolated 12 0
rm(variable, df_variable, lms, pvals, cats, lbls, num_levels, couleurs)
NB : Sous-section qui peut poser problème lors d’un knit si n’est pas executée en même temps que les précédentes (cela ne fonctionne pas bien si les sous-sections précédentes sont en cache).
L’objectif de cette sous-section est d’observer si on observe une différence entre communes fusionnantes et autres communes, du point de vue de l’auto-corrélation spatiale de certains indicateurs.
Globalement, on observe bien des positionnements parfois un peu différents des communes fusionnantes, mais cela est sans doute davantage lié à des effets de contexte (les communes fusionnantes sont sur-représentés dans certains types d’endroits, sous-représentés dans d’autres) que réellement lié à ces clusters.
test_geom$COM_NOUV2 <- dplyr::if_else(test_geom$COM_NOUV == "OUI", 1, 0)
variable <- "COM_NOUV2"
colnames(df2011)
## [1] "CODGEO" "LIBGEO" "CODE_DEPT"
## [4] "CATAEU2010" "REG" "ARR"
## [7] "CV" "UU2010" "AU2010"
## [10] "ZE2010" "EPCI" "P09_ACT1564"
## [13] "P09_CHOM1564" "P09_ETUD1564" "P09_RETR1564"
## [16] "C09_ACT1564_Agr" "C09_ACT1564_ArtCom" "C09_ACT1564_Cadr"
## [19] "C09_ACT1564_ProfInt" "C09_ACT1564_Empl" "C09_ACT1564_Ouvr"
## [22] "P09_EMPLT" "C09_EMPLT_AGRI" "C09_EMPLT_INDUS"
## [25] "C09_EMPLT_CONST" "C09_EMPLT_CTS" "C09_EMPLT_APESAS"
## [28] "P09_ACTOCC" "P09_POP" "P09_POP0014"
## [31] "P09_POP1529" "P09_POP3044" "P09_POP4559"
## [34] "P09_POP6074" "P09_POP75P" "C09_ACTOCC_IN"
## [37] "C09_ACTOCC_OUT" "C09_ACTOCC" "P11_POT_FIN"
## [40] "P11_DGF" "P11_FoyFisc" "P11_Rev_Fisc"
## [43] "P11_IMP_NET" "P11_FoyFisc_Imp" "superficie"
## [46] "CIF_2012" "CIF_2015" "ZAU_POL"
## [49] "ZAU_RUR" "ZAU_MAR_SP" "ZAU_MAR"
## [52] "ZAU_PERI" "ZAU_AU" "FUSION"
## [55] "ChefLieu" "ComDLG" "FusDate"
## [58] "FusPhas" "COM_NOUV" "CODGEO_new"
## [61] "LIBGEO_new" "P09_CHOM1564_RT" "P09_ETUD1564_RT"
## [64] "P09_RETR1564_RT" "C09_ACT1564_Agr_RT" "C09_ACT1564_ArtCom_RT"
## [67] "C09_ACT1564_Cadr_RT" "C09_ACT1564_ProfInt_RT" "C09_ACT1564_Empl_RT"
## [70] "C09_ACT1564_Ouvr_RT" "C09_EMPLT_AGRI_RT" "C09_EMPLT_INDUS_RT"
## [73] "C09_EMPLT_CONST_RT" "C09_EMPLT_CTS_RT" "C09_EMPLT_APESAS_RT"
## [76] "P09_POP0014Y_RT" "P09_POP1529Y_RT" "P09_POP3044Y_RT"
## [79] "P09_POP4559Y_RT" "P09_POP6074Y_RT" "P09_POP75PY_RT"
## [82] "C09_ACTOCC_IN_RT" "C09_ACTOCC_OUT_RT" "C09_EMP_CONC_RT"
## [85] "P11_POT_FIN_RT" "P11_DGF_RT" "P11_Rev_Fisc_RT"
## [88] "P11_IMP_NET_RT" "P11_FoyFisc_Imp_RT"
variable <- "P11_POT_FIN"
variable <- "P09_POP1529Y_RT"
variables_a_tester <- c("P09_CHOM1564", "C09_EMPLT_AGRI","P09_POP0014", "P09_POP1529", "P09_POP6074", "C09_ACTOCC_OUT", "superficie", "P11_POT_FIN")
for (variable in variables_a_tester) {
test_geom2 <- merge(test_geom, df2011[, c("CODGEO", variable)], by = "CODGEO")
colnames(test_geom2)
df_variable <- test_geom2[variable]
lisa <- local_moran(queen_w, df_variable)
# To get the values of the local Moran:
lms <- lisa_values(gda_lisa = lisa)
# To get the pseudo-p values of significance of local Moran computation:
pvals <- lisa_pvalues(lisa)
# To get the cluster indicators of local Moran computation:
cats <- lisa_clusters(lisa, cutoff = 0.05)
# Labels
lbls <- lisa_labels(lisa)
test_geom2$LISA <- lms
test_geom2$pvals <- pvals
test_geom2$cats <- cats
# On rend les catégories plus lisibles en remplaçant le numéro par l'intitulé
test_geom2$cats <- as.factor(test_geom2$cats)
num_levels <- as.numeric(levels(test_geom2$cats))
levels(test_geom2$cats) <- lbls[num_levels + 1]
print(levels(test_geom2$cats))
test_geom2$LISA_retenu <- ifelse(test_geom2$pvals <= 0.1, yes = test_geom2$LISA, no = NA)
par(mar = c(0,0,0,0), mfrow = c(1,2))
couleurs <- c("#f7f7f7", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33")
# Cf. https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6 pour composer les palettes
typoLayer(x = test_geom2 , var = "cats",
# col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
col = couleurs[1:length(levels(test_geom2$cats))],
border = NA,
legend.pos = "topleft",
legend.values.order = levels(test_geom2$cats),
legend.title.txt = paste0("Indices de Moran significatifs\nsur la variable « ", variable, " »\n(pvalue < 0,05)"))
typoLayer(x = test_geom2, var = "COM_NOUV",
legend.values.order = c("OUI", "NON"),
col=c("red","blue"),
border = NA,
legend.pos = "topleft",
legend.title.txt = "Communes fusionnantes")
# plot(test_geom2CN$geometry, col = NA, border = "black", lwd = 1, add = TRUE)
par(mar = c(0,0,0,0), mfrow = c(1,1))
# On regarde la différence entre les communes fusionnantes et les autres
tabcont <- table(test_geom2$cats, test_geom2$COM_NOUV)
print(tabcont) # En valeur absolue
print(round(100*prop.table(tabcont,margin=1),1)) # Pourcentages, le total se fait par lignes
# round(100*prop.table(tabcont,margin=),1) # Pourcentages, le total se fait sur l'ensemble de la population
print(round(100*prop.table(tabcont,margin=2),1)) # Pourcentages, le total se fait par colonnes
test<-chisq.test(tabcont)
# test$observed
# round(test$expected,1)
# round(test$residuals,2)
print(test)
}
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 25521 2148
## High-High 1553 27
## Low-Low 5498 352
## Low-High 937 30
## High-Low 116 14
## Isolated 12 0
##
## NON OUI
## Not significant 92.2 7.8
## High-High 98.3 1.7
## Low-Low 94.0 6.0
## Low-High 96.9 3.1
## High-Low 89.2 10.8
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 75.9 83.5
## High-High 4.6 1.1
## Low-Low 16.3 13.7
## Low-High 2.8 1.2
## High-Low 0.3 0.5
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 125.46, df = 5, p-value < 2.2e-16
##
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 24929 1921
## High-High 2724 294
## Low-Low 4425 235
## Low-High 1188 97
## High-Low 359 24
## Isolated 12 0
##
## NON OUI
## Not significant 92.8 7.2
## High-High 90.3 9.7
## Low-Low 95.0 5.0
## Low-High 92.5 7.5
## High-Low 93.7 6.3
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 74.1 74.7
## High-High 8.1 11.4
## Low-Low 13.2 9.1
## Low-High 3.5 3.8
## High-Low 1.1 0.9
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 63.652, df = 5, p-value = 2.133e-12
##
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 24406 2117
## High-High 1902 38
## Low-Low 6576 374
## Low-High 603 24
## High-Low 138 18
## Isolated 12 0
##
## NON OUI
## Not significant 92.0 8.0
## High-High 98.0 2.0
## Low-Low 94.6 5.4
## Low-High 96.2 3.8
## High-Low 88.5 11.5
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 72.6 82.3
## High-High 5.7 1.5
## Low-Low 19.5 14.5
## Low-High 1.8 0.9
## High-Low 0.4 0.7
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 155.87, df = 5, p-value < 2.2e-16
##
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 24796 2150
## High-High 1661 34
## Low-Low 6315 355
## Low-High 741 21
## High-Low 112 11
## Isolated 12 0
##
## NON OUI
## Not significant 92.0 8.0
## High-High 98.0 2.0
## Low-Low 94.7 5.3
## Low-High 97.2 2.8
## High-Low 91.1 8.9
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 73.7 83.6
## High-High 4.9 1.3
## Low-Low 18.8 13.8
## Low-High 2.2 0.8
## High-Low 0.3 0.4
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 153.54, df = 5, p-value < 2.2e-16
##
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 25292 2180
## High-High 2005 42
## Low-Low 5495 302
## Low-High 676 25
## High-Low 157 22
## Isolated 12 0
##
## NON OUI
## Not significant 92.1 7.9
## High-High 97.9 2.1
## Low-Low 94.8 5.2
## Low-High 96.4 3.6
## High-Low 87.7 12.3
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 75.2 84.8
## High-High 6.0 1.6
## Low-Low 16.3 11.7
## Low-High 2.0 1.0
## High-Low 0.5 0.9
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 161.05, df = 5, p-value < 2.2e-16
##
## [1] "High-High" "Isolated"
##
## NON OUI
## High-High 33625 2571
## Isolated 12 0
##
## NON OUI
## High-High 92.9 7.1
## Isolated 100.0 0.0
##
## NON OUI
## High-High 100 100
## Isolated 0 0
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabcont
## X-squared = 0.15665, df = 1, p-value = 0.6923
##
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 23706 1902
## High-High 3974 279
## Low-Low 4376 264
## Low-High 1228 104
## High-Low 341 22
## Isolated 12 0
##
## NON OUI
## Not significant 92.6 7.4
## High-High 93.4 6.6
## Low-Low 94.3 5.7
## Low-High 92.2 7.8
## High-Low 93.9 6.1
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 70.5 74.0
## High-High 11.8 10.9
## Low-Low 13.0 10.3
## Low-High 3.7 4.0
## High-Low 1.0 0.9
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 22.554, df = 5, p-value = 0.0004106
##
## [1] "Not significant" "High-High" "Low-Low" "Low-High"
## [5] "High-Low" "Isolated"
##
## NON OUI
## Not significant 25369 2177
## High-High 1665 32
## Low-Low 5741 323
## Low-High 749 27
## High-Low 101 12
## Isolated 12 0
##
## NON OUI
## Not significant 92.1 7.9
## High-High 98.1 1.9
## Low-Low 94.7 5.3
## Low-High 96.5 3.5
## High-Low 89.4 10.6
## Isolated 100.0 0.0
##
## NON OUI
## Not significant 75.4 84.7
## High-High 4.9 1.2
## Low-Low 17.1 12.6
## Low-High 2.2 1.1
## High-Low 0.3 0.5
## Isolated 0.0 0.0
##
## Pearson's Chi-squared test
##
## data: tabcont
## X-squared = 144.26, df = 5, p-value < 2.2e-16
rm(variable, df_variable, lms, pvals, cats, lbls, num_levels, couleurs)
Sans lien direct avec l’effet de grappe, éléments de cartographie, qui sont développés dans le chapitre X de la thèse (utilisation du même Markdown pour limiter la création de multiples dépôts).
Cf. l’idée donnée ici : https://www.bnsp.insee.fr/ark:/12148/bc6p08tp83z/f1.pdf#page=2
geom2011 <- st_read("data/geom.gpkg", layer = "geom2011", quiet = TRUE)
geom_new <- st_read("data/geom.gpkg", layer = "geom_new", quiet = TRUE)
geomfus2011 <- st_read("data/geom.gpkg", layer = "geomfus2011", quiet = TRUE)
geomCN_new <- st_read("data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)
dep <- st_read("data/geom.gpkg", layer = "dep", quiet = TRUE)
# Métadonnées
# Liste toutes les variables disponiles
variables_dispo <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ind_target"))
# Liste les variables marquées dans le fichier meta_budgets.xlsx comme nous intéressant
target <- subset(variables_dispo, variable_selec == "X")
Les données socio-économiques qui décrivent les communes en
géographies 2011 et 2021 sont ici importées. On commence par extraire
les communes ayant participé à la création d’une commune nouvelle,
appelées ici communes fusionnantes (datafus2011), les
communes nouvelles, avec les géométries au 1er janvier 2021 et
caractérisées par les données à la géométrie 2011 agrégées
(dataCN_new), ainsi que les communes, à la géométrie 2011,
qui n’ont pas participé à la création d’une commune nouvelle
(dataNfus2011)
load("data/refdata.Rdata")
datafus2011 <- subset(df2011, COM_NOUV == "OUI")
dataCN_new <- subset(df_new, COM_NOUV == "OUI")
dataNfus2011 <- subset(df2011, COM_NOUV == "NON")
Dans un certain nombre de cas, il sera utile d’avoir, dans un même objet, les données et les géométries. Les données sont ici jointes aux couches géographiques d’intérêt.
geom2011 <- merge(geom2011, df2011, by = "CODGEO")
geom_new <- merge(geom_new, df_new, by = "CODGEO_new")
geomCN_new <- merge(geomCN_new, dataCN_new, by = "CODGEO_new")
geomfus2011 <- merge(geomfus2011, datafus2011, by = "CODGEO")
# par(mfrow = c(1,1), mar=c(0,0,1.5,0))
# En fonction de la population de la commune fusionnante
choroLayer(x = geomfus2011 , var = "P09_POP",
method = "quantile", nclass = 4,
col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 2, n2 = 2),
border = NA,
legend.pos = "topleft", legend.values.rnd = 2,
legend.title.txt = "Nombre d'habitants (2009,\nregroupement par quartiles)")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
layoutLayer(title = " ",# "Communes fusionnantes (2011-2024)\nen fonction de leur nombre d'habitants",
author = "G. Bideau, 2024",
tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
sources = "Sources : INSEE, IGN, 2024")
# On crée un tableau donnant, pour chaque département, la population moyenne des communes fusionnantes
tableau_dep <- as.data.frame(tapply(geomfus2011$P09_POP, INDEX = geomfus2011$CODE_DEPT, mean))
colnames(tableau_dep) <- c("moyenne_Cfus")
tableau_dep$median_Cfus <- tapply(geomfus2011$P09_POP, INDEX = geomfus2011$CODE_DEPT, median)
tableau_dep$CODE_DEPT <- row.names(tableau_dep)
tableau_dep$Nbr_Cfus <- table(geomfus2011$CODE_DEP)
tableau_dep <- merge(tableau_dep, dep[, c("CODE_DEPT", "LIBELLE")], by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)
tableau_dep_tt <- as.data.frame(tapply(geom2011$P09_POP, INDEX = geom2011$CODE_DEPT, mean))
colnames(tableau_dep_tt) <- c("moyenne")
tableau_dep_tt$median <- tapply(geom2011$P09_POP, INDEX = geom2011$CODE_DEPT, median)
tableau_dep_tt$CODE_DEPT <- row.names(tableau_dep_tt)
tableau_dep <- merge(tableau_dep, tableau_dep_tt, by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)
tableau_pr_publi <- subset(tableau_dep, tableau_dep$Nbr_Cfus > 40)
tableau_pr_publi <- tableau_pr_publi[order(-tableau_pr_publi$moyenne_Cfus),]
kable(tableau_pr_publi[, c("LIBELLE", "Nbr_Cfus", "moyenne_Cfus", "moyenne", "median_Cfus", "median")],
col.names = c("Département", "Nombre de communes fusionnantes", "Population moyenne des communes fusionnantes", "Population communale moyenne", "Population médiane des communes fusionnantes", "Population communale médiane"), digits = 0)
| Département | Nombre de communes fusionnantes | Population moyenne des communes fusionnantes | Population communale moyenne | Population médiane des communes fusionnantes | Population communale médiane | |
|---|---|---|---|---|---|---|
| 76 | Vendée | 44 | 2455 | 2221 | 1283.0 | 1288.0 |
| 45 | Maine-et-Loire | 225 | 1339 | 2149 | 978.0 | 1030.0 |
| 46 | Manche | 206 | 1128 | 828 | 379.5 | 381.0 |
| 1 | Ain | 45 | 924 | 1405 | 351.0 | 756.0 |
| 68 | Savoie | 48 | 894 | 1348 | 443.5 | 549.0 |
| 73 | Deux-Sèvres | 73 | 862 | 1201 | 466.0 | 579.0 |
| 70 | Seine-Maritime | 48 | 849 | 1678 | 474.5 | 494.0 |
| 24 | Eure | 127 | 720 | 863 | 377.0 | 416.0 |
| 25 | Eure-et-Loir | 54 | 655 | 1056 | 423.5 | 430.0 |
| 12 | Calvados | 221 | 588 | 964 | 317.0 | 340.0 |
| 22 | Doubs | 43 | 587 | 884 | 373.0 | 255.5 |
| 14 | Charente | 64 | 579 | 870 | 356.5 | 387.0 |
| 21 | Dordogne | 78 | 571 | 740 | 270.5 | 359.0 |
| 57 | Orne | 150 | 474 | 579 | 236.0 | 244.0 |
| 80 | Yonne | 43 | 392 | 755 | 307.0 | 354.0 |
| 36 | Jura | 78 | 373 | 480 | 159.0 | 212.5 |
| 44 | Lozère | 50 | 303 | 417 | 179.0 | 181.0 |
# En fonction de la population de la commune nouvelle
choroLayer(x = geomCN_new , var = "P09_POP",
method = "quantile", nclass = 4,
col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 2, n2 = 2),
border = NA,
legend.pos = "topleft", legend.values.rnd = 2,
legend.title.txt = "Nombre d'habitants (2009,\nregroupement par quartiles)")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
plot(st_geometry(geomCN_new), add = TRUE, lwd = 0.2)
layoutLayer(title = "",# "Communes nouvelles (2011-2024)\nen fonction de leur nombre d'habitants",
author = "G. Bideau, 2024",
tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
sources = "Sources : INSEE, IGN, 2024")
# On crée un tableau donnant, pour chaque département, la population moyenne des communes nouvelles
tableau_dep <- as.data.frame(tapply(geomCN_new$P09_POP, INDEX = geomCN_new$CODE_DEPT_new, mean))
colnames(tableau_dep) <- c("moyenne_CN")
tableau_dep$median_CN <- tapply(geomCN_new$P09_POP, INDEX = geomCN_new$CODE_DEPT_new, median)
tableau_dep$CODE_DEPT <- row.names(tableau_dep)
tableau_dep$Nbr_Cfus <- table(geomfus2011$CODE_DEP) # On garde le nombre de communes fusionnantes pour faciliter les comparaisons entre les deux tableaux
tableau_dep <- merge(tableau_dep, dep[, c("CODE_DEPT", "LIBELLE")], by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)
tableau_dep_tt <- as.data.frame(tapply(geom_new$P09_POP, INDEX = geom_new$CODE_DEPT_new, mean))
colnames(tableau_dep_tt) <- c("moyenne")
tableau_dep_tt$median <- tapply(geom_new$P09_POP, INDEX = geom_new$CODE_DEPT_new, median)
tableau_dep_tt$CODE_DEPT <- row.names(tableau_dep_tt)
tableau_dep <- merge(tableau_dep, tableau_dep_tt, by = "CODE_DEPT", all.x = TRUE, all.y = FALSE)
tableau_pr_publi <- subset(tableau_dep, tableau_dep$Nbr_Cfus > 40)
tableau_pr_publi <- tableau_pr_publi[order(-tableau_pr_publi$moyenne_CN),]
kable(tableau_pr_publi[, c("LIBELLE", "Nbr_Cfus", "moyenne_CN", "moyenne", "median_CN", "median")],
col.names = c("Département", "Nombre de communes nouvelles", "Population moyenne des communes nouvelles", "Population communale moyenne", "Population médiane des communes nouvelles", "Population communale médiane"), digits = 0)
| Département | Nombre de communes nouvelles | Population moyenne des communes nouvelles | Population communale moyenne | Population médiane des communes nouvelles | Population communale médiane | |
|---|---|---|---|---|---|---|
| 45 | Maine-et-Loire | 225 | 7922 | 4431 | 5631.5 | 1569.5 |
| 76 | Vendée | 44 | 6279 | 2457 | 3529.0 | 1421.0 |
| 46 | Manche | 206 | 4658 | 1120 | 2508.5 | 444.0 |
| 70 | Seine-Maritime | 48 | 4075 | 1768 | 3228.5 | 503.0 |
| 12 | Calvados | 221 | 3012 | 1289 | 1662.0 | 407.5 |
| 68 | Savoie | 48 | 2684 | 1506 | 2697.5 | 612.0 |
| 73 | Deux-Sèvres | 73 | 2621 | 1431 | 1442.5 | 691.0 |
| 24 | Eure | 127 | 2473 | 996 | 1385.0 | 444.0 |
| 57 | Orne | 150 | 2368 | 759 | 2138.5 | 274.0 |
| 1 | Ain | 45 | 2309 | 1502 | 898.5 | 799.5 |
| 25 | Eure-et-Loir | 54 | 2210 | 1166 | 1522.0 | 472.0 |
| 21 | Dordogne | 78 | 1857 | 819 | 1216.5 | 383.0 |
| 14 | Charente | 64 | 1685 | 971 | 1203.5 | 426.0 |
| 80 | Yonne | 43 | 1533 | 812 | 1067.0 | 381.0 |
| 22 | Doubs | 43 | 1401 | 923 | 1323.5 | 260.0 |
| 36 | Jura | 78 | 1039 | 529 | 825.5 | 235.0 |
| 44 | Lozère | 50 | 890 | 508 | 579.0 | 209.0 |
On regarde la position des communes fusionnantes puis des communes nouvelles vis-à-vis des communes françaises
par(mfrow = c(1,2), mar=c(0,0,1.2,0))
# Création des données concernant les quantiles, mais en s'appuyant sur les données de l'ensemble des communes françaises
geomfus2011$P09_POP_quantiles <- as.factor(cut(geomfus2011$P09_POP,
breaks=c(quantile(geom2011$P09_POP, probs=seq(0, 1, 1/6), na.rm = TRUE))))
geomCN_new$P09_POP_quantiles <- as.factor(cut(geomCN_new$P09_POP,
breaks=c(quantile(geom_new$P09_POP, probs=seq(0, 1, 1/6), na.rm = TRUE))))
# En fonction de la population de la commune fusionnante
typoLayer(x = geomfus2011 , var = "P09_POP_quantiles",
border = NA,
legend.values.order = levels(geomfus2011$P09_POP_quantiles),
col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
# col = carto.pal(pal1 = "red.pal", n1 = 5),
legend.pos = "topleft",
legend.title.txt = "Situation au sein\ndes quantiles des\ncommunes françaises")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
layoutLayer(title = "Communes fusionnantes (2011-2024)\nen fonction de leur nombre d'habitant",
author = "G. Bideau, 2024",
tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
sources = "Sources : INSEE, IGN, 2024")
# En fonction de la population de la commune nouvelle
typoLayer(x = geomCN_new , var = "P09_POP_quantiles",
col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
border = NA,
legend.values.order = levels(geomCN_new$P09_POP_quantiles),
legend.pos = "topleft",
legend.title.txt = "Situation au sein\ndes quantiles des\ncommunes françaises")
plot(st_geometry(dep), add = TRUE, lwd = 0.3)
plot(st_geometry(geomCN_new), add = TRUE, lwd = 0.2)
layoutLayer(title = "Communes nouvelles (2011-2024)\nen fonction de leur nombre d'habitant",
author = "G. Bideau, 2024",
tabtitle = TRUE, frame = FALSE, col = "white", coltitle = "black",
sources = "Sources : INSEE, IGN, 2024")
On voit bien que les communes nouvelles sont bien davantage dans la strate supérieure. Mais il ne faut pas croire non plus que seules de petites communes fusionnent, ce n’est pas le cas.